Evaluation of RISK survival models

This document highlights the use of

for the evaluation (RRPlot), and calibration of cox models (CoxRiskCalibration) or logistic models (CalibrationProbPoissonRisk) of survival data.

Furthermore, it can be used to evaluate any Risk index that reruns the probability of a future event on external data-set.

This document will use the survival::rotterdam, and survival::gbsg data-sets to train and predict the risk of cancer recurrence after surgery. Both Cox and Logistic models will be trained and evaluated.

Here are some sample plots returned by the evaluated functions:

The libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
source("~/GitHub/RISKPLOTS/CODE/auxplots.R")

Breast Cancer Royston-Altman data

data(gbsg, package=“survival”) and data(rotterdam, package=“survival”)

gbsgdata <- gbsg
rownames(gbsgdata) <- gbsgdata$pid
gbsgdata$pid <- NULL

odata <-rotterdam
rownames(odata) <- odata$pid
odata$pid <- NULL
odata$rfstime <- odata$rtime
odata$status <- odata$recur
odata$rtime <- NULL
odata$recur <- NULL

odata <- odata[,colnames(odata) %in% colnames(gbsgdata)]

odata$size <- 10*(odata$size=="<=20") + 
  35*(odata$size=="20-50") + 
  60*(odata$size==">50")

data <- as.data.frame(model.matrix(Surv(rfstime,status)~.*age,odata))

data$`(Intercept)` <- NULL

dataBrestCancerTrain <- cbind(time=odata[rownames(data),"rfstime"],status=odata[rownames(data),"status"],data)

colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),":","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain)," ","")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),"\\.","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),"-","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),">","_")
dataBrestCancerTrain$time <- dataBrestCancerTrain$time/365 ## To years


pander::pander(table(odata[rownames(data),"status"]),caption="rotterdam")
rotterdam
0 1
1464 1518

data(gbsg, package=“survival”) data conditioning

gbsgdata <- gbsgdata[,colnames(odata)]
data <- as.data.frame(model.matrix(Surv(rfstime,status)~.*.,gbsgdata))

data$`(Intercept)` <- NULL

dataBrestCancerTest <- cbind(time=gbsgdata[rownames(data),"rfstime"],status=gbsgdata[rownames(data),"status"],data)

colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),":","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest)," ","")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),"\\.","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),"-","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),">","_")
dataBrestCancerTest$time <- dataBrestCancerTest$time/365

pander::pander(table(odata[rownames(data),"status"]), caption="gbsg")
gbsg
0 1
499 183

Cox Modeling

#ml <- BSWiMS.model(Surv(time,status)~.,data=dataBrestCancerTrain,loops=1,NumberofRepeats = 5)
ml <- BSWiMS.model(Surv(time,status)~.,data=dataBrestCancerTrain,NumberofRepeats = 5)

[++-++-++-++-++-].

sm <- summary(ml)
pander::pander(sm$coefficients)
  Estimate lower HR upper u.Accuracy r.Accuracy full.Accuracy u.AUC r.AUC full.AUC IDI NRI z.IDI z.NRI Delta.AUC Frequency
age_nodes 5.33e-04 1.000 1.001 1.001 0.626 0.588 0.625 0.630 0.588 0.627 0.03069 0.4398 13.21 13.30 0.03894 1
nodes 4.40e-02 1.041 1.045 1.049 0.637 0.599 0.640 0.640 0.600 0.641 0.03197 0.4403 13.10 13.40 0.04147 1
size 7.33e-03 1.005 1.007 1.010 0.595 0.637 0.640 0.595 0.640 0.641 0.01103 0.2963 6.56 8.18 0.00124 1
age_size 7.57e-05 1.000 1.000 1.000 0.567 0.629 0.625 0.568 0.633 0.627 0.00762 0.3044 5.49 8.41 -0.00568 1
grade 2.17e-01 1.148 1.242 1.343 0.565 0.637 0.640 0.561 0.639 0.641 0.00777 0.2420 5.41 7.53 0.00214 1
age -3.91e-03 0.995 0.996 0.998 0.513 0.626 0.640 0.513 0.626 0.641 0.00376 0.0931 4.97 2.54 0.01539 1
age_grade -6.64e-04 0.999 0.999 1.000 0.508 0.616 0.625 0.509 0.619 0.627 0.00186 0.0574 3.44 1.57 0.00843 1
age_pgr -2.60e-06 1.000 1.000 1.000 0.548 0.626 0.625 0.544 0.629 0.627 0.00267 0.1705 2.97 5.21 -0.00183 1

Cox Model Performance

Here we evaluate the model using the RRPlot() function.

The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

timeinterval <- 5 # Five years

h0 <- sum(dataBrestCancerTrain$status & dataBrestCancerTrain$time <= timeinterval)
h0 <- h0/sum((dataBrestCancerTrain$time > timeinterval) | (dataBrestCancerTrain$status==1))

pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
Initial Parameters
h0 timeinterval
0.429 5 
index <- predict(ml,dataBrestCancerTrain)
rdata <- cbind(dataBrestCancerTrain$status,ppoisGzero(index,h0))

rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

By Risk Categories

obsexp <- rrAnalysisTrain$OERatio$atThrEstimates

expObs(obsexp,"Train: Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 1518 1443 1596 1451 1.046 0.994 1.100 8.08e-02
low 819 764 877 879 0.932 0.869 0.998 4.47e-02
90% 225 197 256 188 1.200 1.048 1.367 7.63e-03
80% 474 432 519 382 1.242 1.133 1.359 5.10e-06

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
1738186 5.24e-25 * * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
90068 0.129 two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.364 O:0.364 <= Risk < 0.435 O:High Risk >= 0.435 E:Low Risk < 0.364 E:0.364 <= Risk < 0.435 E:High Risk >= 0.435
1Q 3.42 1.84 1.05 6.34 4.67 2.31
Median 6.83 4.27 2.23 7.55 4.99 3.15
3Q 9.42 7.97 4.78 8.58 5.29 3.86

The Time vs. Events are not calibrated. Lets do the calibration

Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @:0.8 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.435 0.364 0.339 0.239 1.98e-01 0.5000
RR 1.723 1.712 1.777 2.230 9.22e+01 1.7119
RR_LCI 1.617 1.602 1.653 1.839 1.91e-01 1.6049
RR_UCI 1.835 1.829 1.910 2.704 4.46e+04 1.8261
SEN 0.312 0.460 0.586 0.947 1.00e+00 0.2134
SPE 0.899 0.800 0.704 0.172 1.23e-02 0.9426
BACC 0.606 0.630 0.645 0.559 5.06e-01 0.5780
NetBenefit 0.121 0.178 0.224 0.355 3.90e-01 0.0805 
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.05 0.994 1.1 0.0808 
pander::pander(rrAnalysisTrain$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 1518 1443 1596 1451 1.046 0.994 1.100 8.08e-02
low 819 764 877 879 0.932 0.869 0.998 4.47e-02
90% 225 197 256 188 1.200 1.048 1.367 7.63e-03
80% 474 432 519 382 1.242 1.133 1.359 5.10e-06 
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
1.16 1.16 1.15 1.16 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
1.35 1.35 1.35 1.35 
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.678 0.678 0.664 0.692 
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.693 0.674 0.712 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.312 0.289 0.336 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.882 0.914 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.435 0.364 
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 504.157535 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1990 819 1150 95.4 399.4
class=1 370 225 169 18.5 20.9
class=2 622 474 199 382.0 445.9

Cox Calibration

op <- par(no.readonly = TRUE)


calprob <- CoxRiskCalibration(ml,dataBrestCancerTrain,"status","time")

( 7.156166 , 29202.06 , 1.068545 , 1518 , 1795.957 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.714 1.38 7.16

The RRplot() of the calibrated model

rcaldata <- cbind(dataBrestCancerTrain$status,calprob$prob)


rrAnalysisCalTrain <- RRPlot(rcaldata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Cal. Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

By Risk Categories

obsexp <- rrAnalysisCalTrain$OERatio$atThrEstimates

expObs(obsexp,"Cal. Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 1518 1443 1596 1687 0.900 0.855 0.946 2.97e-05
low 819 764 877 1022 0.802 0.748 0.858 6.03e-11
90% 225 197 256 218 1.032 0.902 1.176 6.35e-01
80% 474 432 519 444 1.068 0.974 1.169 1.54e-01

Time to Event Analysis

rrAnalysisdata <- rrAnalysisCalTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
2322576 0.0357 * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
104552 0.0869 two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.529 O:0.529 <= Risk < 0.613 O:High Risk >= 0.613 E:Low Risk < 0.529 E:0.529 <= Risk < 0.613 E:High Risk >= 0.613
1Q 3.42 1.84 1.05 5.46 4.01 1.99
Median 6.83 4.27 2.23 6.49 4.29 2.71
3Q 9.42 7.97 4.78 7.38 4.55 3.32

Calibrated Train Performance

pander::pander(t(rrAnalysisCalTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @:0.8 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.6135 0.529 0.498 0.365 3.07e-01 0.500
RR 1.7227 1.712 1.777 2.230 9.22e+01 1.737
RR_LCI 1.6173 1.602 1.653 1.839 1.91e-01 1.617
RR_UCI 1.8348 1.829 1.910 2.704 4.46e+04 1.866
SEN 0.3123 0.460 0.586 0.947 1.00e+00 0.573
SPE 0.8989 0.800 0.704 0.172 1.23e-02 0.706
BACC 0.6056 0.630 0.645 0.559 5.06e-01 0.640
NetBenefit 0.0802 0.124 0.154 0.249 2.94e-01 0.148 
pander::pander(t(rrAnalysisCalTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
0.9 0.855 0.946 2.97e-05 
pander::pander(t(rrAnalysisCalTrain$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
0.996 0.996 0.989 1 
pander::pander(t(rrAnalysisCalTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
1.05 1.05 1.05 1.05 
pander::pander(rrAnalysisCalTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.678 0.678 0.663 0.692 
pander::pander(t(rrAnalysisCalTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.693 0.674 0.712 
pander::pander((rrAnalysisCalTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.312 0.289 0.336 
pander::pander((rrAnalysisCalTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.882 0.914 
pander::pander(t(rrAnalysisCalTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.613 0.529 
pander::pander(rrAnalysisCalTrain$surdif,caption="Logrank test")
Logrank test Chisq = 504.157535 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1990 819 1150 95.4 399.4
class=1 370 225 169 18.5 20.9
class=2 622 474 199 382.0 445.9

Performance on the external data set

index <- predict(ml,dataBrestCancerTest)
pp <- predictionStats_binary(cbind(dataBrestCancerTest$status,index),plotname="Breast Cancer")

par(op)


prob <- ppoisGzero(index,h0)
rdata <- cbind(dataBrestCancerTest$status,prob)
#rrCoxTestAnalysis <- RRPlot(rdata,atThr=c(0.9,0.8),
rrCoxTestAnalysis <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

External Data Report

pander::pander(t(rrCoxTestAnalysis$keyPoints),caption="Threshold values")
Threshold values
  @:0.435 @:0.364 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.4349 0.364 0.358 0.250 2.12e-01 0.4999
RR 1.7921 1.746 1.773 2.737 2.64e+01 1.7885
RR_LCI 1.5300 1.472 1.490 1.446 5.65e-02 1.5235
RR_UCI 2.0991 2.070 2.110 5.181 1.23e+04 2.0997
SEN 0.3712 0.555 0.585 0.973 1.00e+00 0.2709
SPE 0.8475 0.690 0.667 0.103 1.55e-02 0.9044
BACC 0.6094 0.623 0.626 0.538 5.08e-01 0.5876
NetBenefit 0.0956 0.142 0.150 0.255 2.86e-01 0.0642 
pander::pander(t(rrCoxTestAnalysis$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.38 1.23 1.54 1.36e-07 
pander::pander(rrCoxTestAnalysis$c.index,caption="C. Index")

  • C Index: 0.668

  • Dxy: 0.336

  • S.D.: 0.0309

  • n: 686

  • missing: 0

  • uncensored: 299

  • Relevant Pairs: 266144

  • Concordant: 177775

  • Uncertain: 203702

  • cstatCI:

    mean.C Index median lower upper
    0.668 0.668 0.637 0.698 
pander::pander(t(rrCoxTestAnalysis$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.669 0.629 0.71 
pander::pander((rrCoxTestAnalysis$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.371 0.316 0.429 
pander::pander((rrCoxTestAnalysis$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.85 0.811 0.884 
pander::pander(t(rrCoxTestAnalysis$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.435 0.364 
pander::pander(rrCoxTestAnalysis$surdif,caption="Logrank test")
Logrank test Chisq = 92.126440 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 400 133 198.1 21.376 64.3
class=1 117 55 49.3 0.668 0.8
class=2 169 111 51.7 68.135 84.0

Calibrating the index on the test data

calprob <- CoxRiskCalibration(ml,dataBrestCancerTest,"status","time")

( 5.664115 , 3012.376 , 1.295745 , 299 , 329.4124 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.579 1 5.66 
rdata <- cbind(dataBrestCancerTest$status,calprob$prob)

rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTest$time,
                     title="Cal. Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

By Risk Categories

obsexp <- rrAnalysisTest$OERatio$atThrEstimates

expObs(obsexp,"Cal. Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 299 266.1 334.9 256.6 1.17 1.037 1.30 0.00954
low 171 146.3 198.6 143.7 1.19 1.018 1.38 0.02692
90% 43 31.1 57.9 29.9 1.44 1.041 1.94 0.02174
80% 85 67.9 105.1 80.8 1.05 0.841 1.30 0.61637

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTest

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
26492 3.01e-69 * * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
3703 0.669 two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.508 O:0.508 <= Risk < 0.599 O:High Risk >= 0.599 E:Low Risk < 0.508 E:0.508 <= Risk < 0.599 E:High Risk >= 0.599
1Q 1.95 1.37 1.08 4.95 3.34 1.58
Median 3.35 2.35 1.72 5.72 3.65 2.30
3Q 4.83 3.69 3.01 6.50 3.86 2.69

After Calibration Report

pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @:0.8 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.5985 0.5082 0.450 0.322 2.75e-01 0.5002
RR 1.8121 1.7321 1.773 2.737 2.64e+01 1.7553
RR_LCI 1.5464 1.4746 1.490 1.446 5.65e-02 1.4941
RR_UCI 2.1233 2.0346 2.110 5.181 1.23e+04 2.0620
SEN 0.2876 0.4281 0.585 0.973 1.00e+00 0.4415
SPE 0.8992 0.7959 0.667 0.103 1.55e-02 0.7907
BACC 0.5934 0.6120 0.626 0.538 5.08e-01 0.6161
NetBenefit 0.0406 0.0676 0.101 0.184 2.25e-01 0.0743 
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.17 1.04 1.3 0.00954 
pander::pander(rrAnalysisTest$c.index,caption="C. Index")

  • C Index: 0.668

  • Dxy: 0.336

  • S.D.: 0.0309

  • n: 686

  • missing: 0

  • uncensored: 299

  • Relevant Pairs: 266144

  • Concordant: 177775

  • Uncertain: 203702

  • cstatCI:

    mean.C Index median lower upper
    0.668 0.668 0.638 0.695 
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.669 0.629 0.71 
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.284 0.234 0.339 
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.865 0.927 
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.599 0.508 
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 89.011587 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 479 171 232.5 16.28 74.64
class=1 83 43 30.2 5.38 6.02
class=2 124 85 36.2 65.61 75.76

Logistic Model

Here we train a logistic model on the same data set

## Only label subjects that present event withing five years

dataBrestCancerR <- subset(dataBrestCancerTrain, time>=5 | status==1)
dataBrestCancerR$status <- dataBrestCancerR$status * (dataBrestCancerR$time < 5)
dataBrestCancerR$time <- NULL

#ml <- BSWiMS.model(status~1,data=dataBrestCancerR,loops=20,NumberofRepeats = 5)
#mlog <- BSWiMS.model(status~1,data=dataBrestCancerR,loops=1,NumberofRepeats = 5)
mlog <- BSWiMS.model(status~1,data=dataBrestCancerR,NumberofRepeats = 5)

[++-++-++-++-++-].

sm <- summary(mlog)
pander::pander(sm$coefficients)
  Estimate lower OR upper u.Accuracy r.Accuracy full.Accuracy u.AUC r.AUC full.AUC IDI NRI z.IDI z.NRI Delta.AUC Frequency
age_nodes 0.001336 1.001 1.001 1.001 0.678 0.608 0.684 0.642 0.582 0.653 0.07150 0.5431 14.27 15.33 0.07040 1
nodes 0.086277 1.081 1.090 1.099 0.676 0.632 0.691 0.639 0.622 0.663 0.07210 0.5673 14.21 16.08 0.04085 1
size 0.010294 1.007 1.010 1.014 0.611 0.686 0.691 0.618 0.653 0.663 0.01523 0.3186 6.26 8.38 0.00992 1
age_size 0.000151 1.000 1.000 1.000 0.608 0.678 0.684 0.577 0.643 0.653 0.01434 0.3004 6.09 7.89 0.00993 1
grade 0.375662 1.280 1.456 1.656 0.571 0.683 0.691 0.500 0.653 0.663 0.01179 0.1809 5.72 4.76 0.01037 1
age_meno -0.003322 0.995 0.997 0.998 0.571 0.684 0.684 0.500 0.652 0.653 0.00828 0.0799 4.77 2.08 0.00126 1
pgr -0.000254 1.000 1.000 1.000 0.571 0.681 0.684 0.500 0.649 0.653 0.00398 0.1921 3.31 5.57 0.00381 1
age_grade -0.001866 0.997 0.998 1.000 0.574 0.690 0.691 0.507 0.662 0.663 0.00273 0.0713 2.62 1.85 0.00127 1

Logistic Model Performance

op <- par(no.readonly = TRUE)


cprob <- predict(mlog,dataBrestCancerTrain)

rdata <- cbind(dataBrestCancerTrain$status,cprob)
rrAnalysisLogTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Logistic Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=5.0)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogTrain$OERatio$atThrEstimates

expObs(obsexp,"Logistic: Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 1518 1443 1596 1857 0.817 0.777 0.860 4.90e-16
low 814 759 872 1006 0.809 0.754 0.866 4.24e-10
90% 215 187 246 253 0.850 0.740 0.972 1.68e-02
80% 489 447 534 594 0.824 0.752 0.900 1.11e-05

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
2352145 0.00635 * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
136255 4.55e-14 * * * two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Logistic: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.451 O:0.451 <= Risk < 0.560 O:High Risk >= 0.560 E:Low Risk < 0.451 E:0.451 <= Risk < 0.560 E:High Risk >= 0.560
1Q 3.38 2.02 1.06 5.36 3.37 1.41
Median 6.87 4.49 2.23 6.49 3.67 2.01
3Q 9.44 7.97 4.70 8.20 4.04 2.54

Training Report

pander::pander(t(rrAnalysisLogTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @:0.8 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.560 0.451 0.379 0.240 1.88e-01 0.500
RR 1.753 1.719 1.816 2.163 9.73e+01 1.754
RR_LCI 1.647 1.609 1.684 1.725 2.01e-01 1.646
RR_UCI 1.866 1.837 1.958 2.713 4.71e+04 1.868
SEN 0.322 0.464 0.627 0.962 1.00e+00 0.379
SPE 0.900 0.799 0.671 0.122 1.30e-02 0.867
BACC 0.611 0.631 0.649 0.542 5.06e-01 0.623
NetBenefit 0.101 0.155 0.221 0.353 3.97e-01 0.128 
pander::pander(t(rrAnalysisLogTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
0.817 0.777 0.86 4.9e-16 
pander::pander(rrAnalysisLogTrain$c.index,caption="C. Index")

  • C Index: 0.678

  • Dxy: 0.356

  • S.D.: 0.014

  • n: 2982

  • missing: 0

  • uncensored: 1518

  • Relevant Pairs: 6184528

  • Concordant: 4193543

  • Uncertain: 2703838

  • cstatCI:

    mean.C Index median lower upper
    0.678 0.679 0.665 0.692 
pander::pander(t(rrAnalysisLogTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.694 0.676 0.713 
pander::pander((rrAnalysisLogTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.322 0.299 0.346 
pander::pander((rrAnalysisLogTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.9 0.883 0.915 
pander::pander(t(rrAnalysisLogTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.56 0.451 
pander::pander(rrAnalysisLogTrain$surdif,caption="Logrank test")
Logrank test Chisq = 523.760672 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1984 814 1146 96.4 399.4
class=1 362 215 169 12.4 13.9
class=2 636 489 202 405.9 475.6

Results on the validation set using Logistic model

pre <- predict(mlog,dataBrestCancerTest)
rdata <- cbind(dataBrestCancerTest$status,pre)

rrAnalysisLogTest <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Logistic Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=5)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogTest$OERatio$atThrEstimates

expObs(obsexp,"Logistic Test: Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 299 266.1 334.9 293.4 1.019 0.907 1.14 0.726
low 71 55.5 89.6 59.9 1.186 0.926 1.50 0.155
90% 46 33.7 61.4 44.8 1.027 0.752 1.37 0.822
80% 182 156.5 210.4 186.2 0.978 0.841 1.13 0.798

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogTest

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
44645 4.21e-45 * * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
24672 0.327 two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Logistic Test: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.364 O:0.364 <= Risk < 0.435 O:High Risk >= 0.435 E:Low Risk < 0.364 E:0.364 <= Risk < 0.435 E:High Risk >= 0.435
1Q 2.23 1.70 1.30 6.07 4.63 1.83
Median 3.70 3.24 2.27 6.61 4.92 2.70
3Q 4.98 4.72 3.77 7.55 5.28 3.69

Validation Report

pander::pander(t(rrAnalysisLogTest$keyPoints),caption="Threshold values")
Threshold values
  @:0.435 @:0.364 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.435 0.364 0.452 0.289 2.17e-01 0.5032
RR 1.728 1.578 1.771 2.500 1.75e+01 1.7144
RR_LCI 1.447 1.274 1.491 1.520 3.85e-02 1.4554
RR_UCI 2.063 1.953 2.103 4.113 7.98e+03 2.0195
SEN 0.609 0.763 0.572 0.957 1.00e+00 0.4783
SPE 0.630 0.401 0.680 0.147 1.03e-02 0.7519
BACC 0.620 0.582 0.626 0.552 5.05e-01 0.6151
NetBenefit 0.105 0.139 0.100 0.221 2.82e-01 0.0668 
pander::pander(t(rrAnalysisLogTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.02 0.907 1.14 0.726 
pander::pander(rrAnalysisLogTest$c.index,caption="C. Index")

  • C Index: 0.663

  • Dxy: 0.327

  • S.D.: 0.0311

  • n: 686

  • missing: 0

  • uncensored: 299

  • Relevant Pairs: 266144

  • Concordant: 176569

  • Uncertain: 203702

  • cstatCI:

    mean.C Index median lower upper
    0.663 0.665 0.634 0.695 
pander::pander(t(rrAnalysisLogTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.664 0.624 0.705 
pander::pander((rrAnalysisLogTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.609 0.551 0.664 
pander::pander((rrAnalysisLogTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.633 0.583 0.681 
pander::pander(t(rrAnalysisLogTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.435 0.364 
pander::pander(rrAnalysisLogTest$surdif,caption="Logrank test")
Logrank test Chisq = 55.976567 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 227 71 117.2 18.19 30.16
class=1 135 46 62.4 4.33 5.48
class=2 324 182 119.4 32.81 55.25

Logistic Model Poisson Calibration

riskdata <- cbind(dataBrestCancerTrain$status,predict(mlog,dataBrestCancerTrain,type="prob"),dataBrestCancerTrain$time)
calprob <- CalibrationProbPoissonRisk(riskdata)

( 7.467216 , 29202.06 , 1.11499 , 1518 , 1835.735 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Logistic Calibration Parameters")
h0 Gain DeltaTime
0.688 1.33 7.47 
timeinterval <- calprob$timeInterval;
gain <- calprob$hazardGain

rdata <- cbind(dataBrestCancerTrain$status,calprob$prob)


rrAnalysisLogCalTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Cal. Logistic Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogCalTrain$OERatio$atThrEstimates

expObs(obsexp,"Logistic Cal: Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 1518 1443 1596 1653 0.918 0.873 0.966 0.000821
low 814 759 872 896 0.909 0.847 0.973 0.006134
90% 215 187 246 225 0.955 0.832 1.092 0.526534
80% 489 447 534 528 0.926 0.845 1.012 0.089670

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogCalTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
1928010 3.15e-10 * * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
125277 2.27e-07 * * * two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Logistic Cal: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.549 O:0.549 <= Risk < 0.664 O:High Risk >= 0.664 E:Low Risk < 0.549 E:0.549 <= Risk < 0.664 E:High Risk >= 0.664
1Q 3.38 2.02 1.06 6.02 3.78 1.59
Median 6.87 4.49 2.23 7.29 4.12 2.26
3Q 9.44 7.97 4.70 9.21 4.54 2.85

Report of the calibrated logistic: training

pander::pander(t(rrAnalysisLogCalTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @:0.8 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.6646 0.549 0.469 0.306 2.42e-01 0.500
RR 1.7529 1.719 1.816 2.163 9.73e+01 1.768
RR_LCI 1.6467 1.609 1.684 1.725 2.01e-01 1.645
RR_UCI 1.8661 1.837 1.958 2.713 4.71e+04 1.900
SEN 0.3221 0.464 0.627 0.962 1.00e+00 0.583
SPE 0.8996 0.799 0.671 0.122 1.30e-02 0.705
BACC 0.6109 0.631 0.649 0.542 5.06e-01 0.644
NetBenefit 0.0663 0.116 0.177 0.300 3.55e-01 0.152 
pander::pander(t(rrAnalysisLogCalTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
0.918 0.873 0.966 0.000821 
pander::pander(rrAnalysisLogCalTrain$c.index,caption="C. Index")

  • C Index: 0.678

  • Dxy: 0.356

  • S.D.: 0.014

  • n: 2982

  • missing: 0

  • uncensored: 1518

  • Relevant Pairs: 6184528

  • Concordant: 4193541

  • Uncertain: 2703838

  • cstatCI:

    mean.C Index median lower upper
    0.678 0.678 0.664 0.692 
pander::pander(t(rrAnalysisLogCalTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.694 0.676 0.713 
pander::pander((rrAnalysisLogCalTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.322 0.299 0.346 
pander::pander((rrAnalysisLogCalTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.9 0.883 0.915 
pander::pander(t(rrAnalysisLogCalTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.664 0.549 
pander::pander(rrAnalysisLogCalTrain$surdif,caption="Logrank test")
Logrank test Chisq = 523.760672 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1984 814 1146 96.4 399.4
class=1 362 215 169 12.4 13.9
class=2 636 489 202 405.9 475.6 
probLog <- predict(mlog,dataBrestCancerTest)
aprob <- adjustProb(probLog,gain)

rdata <- cbind(dataBrestCancerTest$status,aprob)
rrAnalysisTestLogistic <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Cal. Logistic Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

By Risk Categories

obsexp <- rrAnalysisTestLogistic$OERatio$atThrEstimates

expObs(obsexp,"Logistic Test: Expected vs. Observed")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 299 266.07 334.9 261.1 1.145 1.019 1.28 0.0203
low 13 6.92 22.2 14.6 0.892 0.475 1.52 0.7934
90% 45 32.82 60.2 30.6 1.469 1.071 1.97 0.0142
80% 241 211.53 273.4 213.6 1.128 0.990 1.28 0.0645

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTestLogistic

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime and rrAnalysisdata$timetoEventData$expectedTime
Test statistic P value Alternative hypothesis
29769 1.69e-64 * * * two.sided 
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction: rrAnalysisdata$timetoEventData$eTime[highrisk] and rrAnalysisdata$timetoEventData$expectedTime[highrisk]
Test statistic P value Alternative hypothesis
27539 2.53e-25 * * * two.sided 
timesdata <- expObsTime(rrAnalysisdata,title="Logistic Test: Expected vs Observed")

pander::pander(timesdata)
  O:Low Risk < 0.364 O:0.364 <= Risk < 0.435 O:High Risk >= 0.435 E:Low Risk < 0.364 E:0.364 <= Risk < 0.435 E:High Risk >= 0.435
1Q 2.54 2.19 1.45 8.68 6.87 2.57
Median 4.02 3.67 2.51 9.33 7.13 4.16
3Q 5.48 4.66 4.38 10.07 7.70 5.31

Report of the calibrated validation

pander::pander(t(rrAnalysisTestLogistic$keyPoints),caption="Threshold values")
Threshold values
  @:0.435 @:0.364 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.4350 0.364 0.5503 0.365 2.77e-01 0.4997
RR 1.6861 2.421 1.7707 2.500 1.75e+01 1.7332
RR_LCI 1.3315 1.474 1.4907 1.520 3.85e-02 1.4324
RR_UCI 2.1352 3.976 2.1033 4.113 7.98e+03 2.0970
SEN 0.8094 0.957 0.5719 0.957 1.00e+00 0.6789
SPE 0.3566 0.142 0.6796 0.147 1.03e-02 0.5504
BACC 0.5830 0.549 0.6257 0.552 5.05e-01 0.6147
NetBenefit 0.0733 0.140 0.0282 0.141 2.22e-01 0.0426 
pander::pander(t(rrAnalysisTestLogistic$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.15 1.02 1.28 0.0203 
pander::pander(rrAnalysisTestLogistic$c.index,caption="C. Index")

  • C Index: 0.663

  • Dxy: 0.327

  • S.D.: 0.0311

  • n: 686

  • missing: 0

  • uncensored: 299

  • Relevant Pairs: 266144

  • Concordant: 176569

  • Uncertain: 203702

  • cstatCI:

    mean.C Index median lower upper
    0.663 0.664 0.635 0.693 
pander::pander(t(rrAnalysisTestLogistic$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.664 0.624 0.705 
pander::pander((rrAnalysisTestLogistic$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.806 0.757 0.849 
pander::pander((rrAnalysisTestLogistic$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.357 0.309 0.407 
pander::pander(t(rrAnalysisTestLogistic$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.435 0.364 
pander::pander(rrAnalysisTestLogistic$surdif,caption="Logrank test")
Logrank test Chisq = 33.341167 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 68 13 39.2 17.47 20.26
class=1 128 45 63.5 5.39 6.86
class=2 490 241 196.3 10.16 29.83

Comparing the COX and Logistic Models on the Independent Data

pander::pander(t(rrCoxTestAnalysis$OAcum95ci))
mean 50% 2.5% 97.5%
1.24 1.24 1.24 1.25 
pander::pander(t(rrAnalysisTestLogistic$OAcum95ci))
mean 50% 2.5% 97.5%
0.928 0.928 0.925 0.931 
pander::pander(t(rrCoxTestAnalysis$OE95ci))
mean 50% 2.5% 97.5%
1.17 1.17 1.14 1.2 
pander::pander(t(rrAnalysisTestLogistic$OE95ci))
mean 50% 2.5% 97.5%
0.972 0.973 0.948 0.996 
maxobs <- sum(dataBrestCancerTest$status)

par(mfrow=c(1,2),cex=0.75)

plot(rrCoxTestAnalysis$CumulativeOvs[,1:2],type="l",lty=1,
     main="Cumulative Probability",
     xlab="Observed",
     ylab="Cumulative Probability",
     ylim=c(0,maxobs),
     xlim=c(0,maxobs))
lines(rrAnalysisTestLogistic$CumulativeOvs[,1:2],lty=2,col="red")
lines(x=c(0,maxobs),y=c(0,maxobs),lty=3,col="gray")
legend("topleft",legend = c("Cox","Logistic","Ideal"),
       col=c("black","red","gray"),
       lty=c(1,2,3),
       cex=0.75
)

dxcox <- rrCoxTestAnalysis$CumulativeOvs$Cumulative-
       rrCoxTestAnalysis$CumulativeOvs$Observed

dxlogit <- rrAnalysisTestLogistic$CumulativeOvs$Cumulative-
       rrAnalysisTestLogistic$CumulativeOvs$Observed

miny <- min(c(dxcox,dxlogit))
maxy <- max(c(dxcox,dxlogit))
plot(rrCoxTestAnalysis$CumulativeOvs$Observed,
     dxcox,
     main="Cumulative Risk Difference",
     xlab="Observed",
     ylab="Cumulative Risk - Observed",
     type="l",
     ylim=c(miny,maxy),
     lty=1)
lines(rrAnalysisTestLogistic$CumulativeOvs$Observed,
     dxlogit,
     lty=2,
     col="red")
legend("topleft",legend = c("Cox","Logistic"),
       col=c("black","red"),
       lty=c(1,2),
       cex=0.75
)

plot(rrCoxTestAnalysis$OEData[,2:3],type="l",lty=1,
     main="Expected over Time",
     xlab="Observed",
     ylab="Expected",
     ylim=c(0,maxobs),
     xlim=c(0,maxobs))
lines(rrAnalysisTestLogistic$OEData[,2:3],lty=2,col="red")
lines(x=c(0,maxobs),y=c(0,maxobs),lty=3,col="gray")
legend("topleft",legend = c("Cox","Logistic","Ideal"),
       col=c("black","red","gray"),
       lty=c(1,2,3),
       cex=0.75
)

coxdif <- rrCoxTestAnalysis$OEData$Expected-
       rrCoxTestAnalysis$OEData$Observed

logdif <- rrAnalysisTestLogistic$OEData$Expected-
       rrAnalysisTestLogistic$OEData$Observed

miny <- min(c(coxdif,logdif))
maxy <- max(c(coxdif,logdif))

plot(rrCoxTestAnalysis$OEData$Observed,
     coxdif,
      ylim=c(miny,maxy),

     main="Expected vs Observed Difference",
     xlab="Observed",
     ylab="Cumulative - Observed",
     type="l",
     lty=1)
lines(rrAnalysisTestLogistic$OEData$Observed,
     logdif,
     lty=2,col="red")

legend("bottomleft",legend = c("Cox","Logistic"),
       col=c("black","red"),
       lty=c(1,2),
       cex=0.75
)

par(op)